%-------------------------------------------------------------------------;
% Horizon Bias and the Term Structure of Equity Returns;
% This code replicates Tables 4 and 7;
%-------------------------------------------------------------------------;

clear all;
clc;
close all;

%-----------------------------;
% Import term structure data;
%-----------------------------;

filename='Term_structure_variables';
data_ts=xlsread(filename);

b=size(data_ts,1); %Length of the data;

Div_ret=data_ts(2:end,6); %BBK extended;
SP_ret=data_ts(2:end,4);

Div_ret=log(1+Div_ret);
SP_ret=log(1+SP_ret);
                                                 
etp=SP_ret-Div_ret;

SP_price=data_ts(1:end,1);
SP_Div=data_ts(1:end,3);
SP_strip_18=data_ts(1:end,5);

SP_pd=log(SP_price./SP_Div);
SP_pd_strip=log(SP_strip_18./SP_Div);

SP_pd=SP_pd(1:b-12);
SP_pd_strip=SP_pd_strip(1:b-12);

%-------------------------------;
% Import explanatory variables;
%-------------------------------;

filename='Explanatory_variables';
data_ev=xlsread(filename);

data_ev=data_ev(1:b-12,:);

year=data_ev(1:end,1);
month=data_ev(1:end,2);

Ltg=data_ev(1:end,3);   
Stg=data_ev(1:end,4);   

%Take logs;
Ltg=log(1+Ltg);
Stg=log(1+Stg);

%Unaffiliated;
Ltg_un=data_ev(1:end,5);
Stg_un=data_ev(1:end,6); 

Ltg_un=log(1+Ltg_un);
Stg_un=log(1+Stg_un);

%Skewness;
skew1=data_ev(1:end,7);
skew1=(skew1-mean(skew1))./std(skew1);

skew5=data_ev(1:end,8);
skew5=(skew5-mean(skew5))./std(skew5);
Diff_skew=skew5-skew1;

%Uncertainty;
uncert1=data_ev(1:end,9);
uncert1=(uncert1-mean(uncert1))./std(uncert1);

uncert5=data_ev(1:end,10);
uncert5=(uncert5-mean(uncert5))./std(uncert5);
Diff_uncert=uncert5-uncert1;

%Other variables capturing business cycle;
rec=data_ev(1:end,11);
default=data_ev(1:end,12);
term=data_ev(1:end,13);
cay=data_ev(1:end,14);

%Horizon bias;
acorr=autocorr(Stg);
theta_0=mean(Stg);
theta_1=acorr(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg;

Stg_impl=Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two;

HB=Ltg-Stg_impl;

%Horizon bias - unaffiliated;

acorr=autocorr(Stg_un);
theta_0=mean(Stg_un);
theta_1=acorr(13,1);

Stg_impl_one=(1/5)*(4*theta_0+3*theta_0*theta_1+2*theta_0*(theta_1^2)+theta_0*(theta_1^3));
Stg_impl_two=(1/5)*(1+theta_1+theta_1^2+theta_1^3+theta_1^4)*Stg_un;

Stg_impl_un=Stg_impl_one*ones(size(Stg_impl_one,1),1)+Stg_impl_two;

HB_un=Ltg_un-Stg_impl_un;

%-------------------------;
% Predictive regressions;
%-------------------------;

h=12;  %Set horizon equal to 12 months;

%Number of lags for overlapping observations;
lags=h-1;

%Sum returns over h periods;
ret_SP=filter(ones(h,1),1,SP_ret);
ret_Div=filter(ones(h,1),1,Div_ret);

ret_SP=ret_SP(h:end);
ret_Div=ret_Div(h:end);

Ltg_h=Ltg(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
Stg_h=Stg(1:min(size(Stg,1),(size(Stg,1)-h+12)));

Ltg_h_un=Ltg_un(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
Stg_h_un=Stg_un(1:min(size(Stg,1),(size(Stg,1)-h+12)));

ret_SP=ret_SP(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
ret_Div=ret_Div(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

SP_pd=SP_pd(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
SP_pd_strip=SP_pd_strip(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

HB=HB(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
HB_un=HB_un(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

rec=rec(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
default=default(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
term=term(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
cay=cay(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

skew1=skew1(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));
skew5=skew5(1:min(size(Ltg,1),(size(Ltg,1)-h+12)));

%Dependent variable;
y=(ret_SP-ret_Div);

%------------------------------------------------;
%Table 4;
%------------------------------------------------;

%Column 1: HB;

x=[ones(size(Ltg_h,1),1) HB];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table4(1:2,1)=[b(2,1);tb(2,1)];
Table4(15:16,1)=[length(error);AR2];

%Column 2: HB + controls;

x=[ones(size(Ltg_h,1),1) HB SP_pd term default cay];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table4(1:10,2)=[b(2,1);tb(2,1);b(3,1);tb(3,1);b(4,1);tb(4,1);b(5,1);tb(5,1);b(6,1);tb(6,1)];
Table4(15:16,2)=[length(error);AR2];

%Column 3: HB + controls + error;

x=[ones(size(Ltg_h,1),1) HB];
[b,b_std,tb,R2,AR2,error]=olsnw((SP_pd-SP_pd_strip),x,1,lags);

x=[ones(size(Ltg_h,1),1) HB SP_pd term default cay error];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table4(1:12,3)=[b(2,1);tb(2,1);b(3,1);tb(3,1);b(4,1);tb(4,1);b(5,1);tb(5,1);b(6,1);tb(6,1);b(7,1);tb(7,1)];
Table4(15:16,3)=[length(error);AR2];

%Column 4: HB + controls + error + NBER;

x=[ones(size(Ltg_h,1),1) HB];
[b,b_std,tb,R2,AR2,error]=olsnw((SP_pd-SP_pd_strip),x,1,lags);

x=[ones(size(Ltg_h,1),1) HB SP_pd term default cay error rec];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table4(1:14,4)=[b(2,1);tb(2,1);b(3,1);tb(3,1);b(4,1);tb(4,1);b(5,1);tb(5,1);b(6,1);tb(6,1);b(7,1);tb(7,1);b(8,1);tb(8,1)];
Table4(15:16,4)=[length(error);AR2];

%Column 5: HB_un;

x=[ones(size(Ltg_h,1),1) HB_un];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table4(1:2,5)=[b(2,1);tb(2,1)];
Table4(15:16,5)=[length(error);AR2];

%Column 6: HB_un + controls;

x=[ones(size(Ltg_h,1),1) HB_un SP_pd term default cay];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table4(1:10,6)=[b(2,1);tb(2,1);b(3,1);tb(3,1);b(4,1);tb(4,1);b(5,1);tb(5,1);b(6,1);tb(6,1)];
Table4(15:16,6)=[length(error);AR2];

%Column 7: HB_un + controls + error;

x=[ones(size(Ltg_h,1),1) HB_un];
[b,b_std,tb,R2,AR2,error]=olsnw((SP_pd-SP_pd_strip),x,1,lags);

x=[ones(size(Ltg_h,1),1) HB_un SP_pd term default cay error];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table4(1:12,7)=[b(2,1);tb(2,1);b(3,1);tb(3,1);b(4,1);tb(4,1);b(5,1);tb(5,1);b(6,1);tb(6,1);b(7,1);tb(7,1)];
Table4(15:16,7)=[length(error);AR2];

%Column 8: HB_un + controls + error + NBER;

x=[ones(size(Ltg_h,1),1) HB_un];
[b,b_std,tb,R2,AR2,error]=olsnw((SP_pd-SP_pd_strip),x,1,lags);

x=[ones(size(Ltg_h,1),1) HB_un SP_pd term default cay error rec];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table4(1:14,8)=[b(2,1);tb(2,1);b(3,1);tb(3,1);b(4,1);tb(4,1);b(5,1);tb(5,1);b(6,1);tb(6,1);b(7,1);tb(7,1);b(8,1);tb(8,1)];
Table4(15:16,8)=[length(error);AR2]

%------------------------------------------------;
%Table 7;
%------------------------------------------------;

%Column 1: HB;

x=[ones(size(Ltg_h,1),1) HB];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table7(1:2,1)=[b(2,1);tb(2,1)];
Table7(11:12,1)=[length(error);AR2];

%Column 2: Skew5-Skew1;

x=[ones(size(Ltg_h,1),1) skew5-skew1];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table7(3:4,2)=[b(2,1);tb(2,1)];
Table7(11:12,2)=[length(error);AR2];

%Column 3: HB + Skew5-Skew1;

x=[ones(size(Ltg_h,1),1) HB (skew5-skew1) HB.*(skew5-skew1)];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table7(1:6,3)=[b(2,1);tb(2,1);b(3,1);tb(3,1);b(4,1);tb(4,1)];
Table7(11:12,3)=[length(error);AR2];

%Column 5: uncert5-uncert1;

x=[ones(size(Ltg_h,1),1) uncert5-uncert1];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table7(7:8,4)=[b(2,1);tb(2,1)];
Table7(11:12,4)=[length(error);AR2];

%Column 6: HB + uncert5-uncert1;

x=[ones(size(Ltg_h,1),1) HB (uncert5-uncert1) HB.*(uncert5-uncert1)];
[b,b_std,tb,R2,AR2,error]=olsnw(y,x,1,lags);

Table7(1:2,5)=[b(2,1);tb(2,1)];
Table7(7:10,5)=[b(3,1);tb(3,1);b(4,1);tb(4,1)];
Table7(11:12,5)=[length(error);AR2]


